#################################################################################
#                                                                               
#                        MEAN CUMULATIVE COUNT                                  
#                                                                               
# Last updated July, 2025.  
#################################################################################
###### NOTE:   The code replies on the function "crprep" which came from the  "mstate" package but it worked perfect on R 3.5.X or eralier versions. With later R version, it runs without any error but the results are not as expected. So please use R 3.5.X version from https://cran.r-project.org/bin/windows/base/old/   
#################################################################################                                           
#################################################################################
# R code to accompany:
#
#   Dong H, Robison LL, Leisenring WM, Martin LJ, Armstrong GT, Yasui Y. (in press).
#   Estimating the burden of recurrent events in the presence of competing risks: 
#   The method of mean cumulative count. American Journal of Epidemiology. 
#
# NOTE:--------------------------------------------------------------------------#
#	This code is to calculate MCC through the Sum of Cumualtive Incidence (SCI), with the option for 95% confidence interval.
# -------------------------------------------------------------------------------#
#
# R version 3.2.4

# Contact: Qi Liu, ql3@ualberta.ca

############ -------------  PARAMETERS -------------- #########
############----- id:    Participants's ID
############----- time:  Follow up time of participant to interest/competing-risk event/censoring
############----- cause: Event of interest is coded as 1; competing-risk event as 2; censoring as 0 
############----- Tstart: Vector with start time of follow-up. Default is 0 (no left truncation). 
############----- type:  Character string specifying the type of calcultion. 
#				 Possible values are "MCC" (MCC calculation by equation) and "SCI" (sum of cumulative incidence). 


#################################################################################


library("mstate")  ###cumulative incidence that can be used for both left truncation and/or right censoring
library(cmprsk)  ###cumulative incidence for right censoring only
### These two packages are needed to fill in NA with previous non-missing value ###
library(dplyr)
#library(DataCombine) #FillDown: Fills in missing (NA) values with the previous non-missing value from "DataCombine" package. not in R anu more.
library(zoo)


###  Revised version of function "crprep" from "mstate" package--Found a bug and contact the author of "mstate". The author sent me the revised version of this function and the package will be updated by the author later.
source("R:\\Biostatistics\\Biostatistics\\Common\\Qi\\webpage\\MCC\\crprep_Revised.R")

###################   This is the main function to run  ###################

MCC=function(id,time,cause,Tstart=0,type="SCI",c_interval=FALSE)
{
	if(type=="MCC"){
		outdata=MCC_eq(id,time,cause)
	}
	if(type=="SCI" & !c_interval){
		outdata=SCI(id,time,cause,Tstart)
	}
	if(type=="SCI" & c_interval){
		outdata=SCI_95CI(id,time,cause,Tstart)
	}
return(outdata)
}



###################  Calculate the Sum of Cumulative Incidence  ########## 

SCI=function(id,time,cause,Tstart)
{
### if the number of event 1 is 0, then no need to run, but give MCC =0;
if(sum(cause==1)==0){
	MCC.final=matrix(c(20,0),nrow=1)
	colnames(MCC.final)=c("time","SumCIs")
}

if(sum(cause==1)>0){
if(max(Tstart)==0) Tstart=rep(0,length(time)) ### By default it is a value 0, need to make it a vector of the same lengnth
### if event time is at cohrot entry, then add a small number to make them different, since in coutning process (to,t1] where t1>t0.
time[time==Tstart]=Tstart[time==Tstart]+exp(-13) 

#### one person ends with either 2 (competing risk) or 0 (censoring), cannot have both 2 and 0.
### if one person had an event and then censoring or competing risk on the same date, put event first.
indata=data.frame(id,time,cause,Tstart)
indata$cause1=indata$cause
indata$cause1[indata$cause==2]=-1
ii=order(indata$id,indata$time,-indata$cause1) 
indata=indata[ii,]

  ## Check whether Tstart is needed. This indicator will be used in SCI fucntion.
  calc.trunc <- any(Tstart[!is.na(Tstart)] != 0)

indata$first <- ave(indata$time, indata$id, FUN = seq_along)
M=max(indata$first) 		#maximum numer of events

### Take the last row so we know the maximum number of events per person
event.number=do.call(rbind, lapply(split(indata, indata$id), tail, 1))[,c("id","first")]
colnames(event.number)=c("id","maxE")

alldata=merge(indata,event.number,by.x="id",by.y="id")

#### make data for cumulative incidence, with dimension M*N
data.MCC=NULL
for(i in 1:M){
	if(i==1){
		data_temp=alldata[alldata$first==1,]
		data_temp$m_event=i
		data.MCC=rbind(data.MCC,data_temp)
		}
	if(i>1){
		## for those with i or more records, take the ith record
		the_ith=alldata[alldata$first==i,]
		##for those with <i records, take the last record. If the last record is an event 1, then change it to 0;
		## because it contributed in data (i-1) already.
		the_last=alldata[alldata$maxE<i & alldata$first==alldata$maxE,]
		the_last$cause[the_last$cause==1]=0
		CIdata_ith=rbind(the_ith,the_last)
		CIdata_ith$m_event=i
		data.MCC=rbind(data.MCC,CIdata_ith)
		}
}

#Calculate cumulative incidence for M times event#
#Store the results in matrix MCC.base#
MCC.base=NULL
for (j in 1:M)
{
	dataj=data.MCC[data.MCC$m_event==j,]
	if(1 %in% dataj$cause){  ##only if the data contains the event of interest
		if (calc.trunc) {
			expdata=crprep(Tstop="time", status="cause", data=dataj, trans = 1, cens = 0, Tstart="Tstart", id="id")
			fit1=survfit(Surv(Tstart,Tstop,status==1)~1,data=expdata,weight=weight.cens*weight.trunc)
		 	cm1=data.frame(Time=summary(fit1)$time,cm=1-summary(fit1)$surv)
		} else { 
			ciGrey=cuminc(dataj$time,dataj$cause) 
			cm1=data.frame(Time=ciGrey[[1]]$time,cm=ciGrey[[1]]$est)
		}
		 cm1$Deta=c(cm1$cm[1],diff(cm1$cm))
		 cm1$cumI=j
		MCC.base=rbind(MCC.base,cm1)
	}
}

#Only keep the time points that affect MCC#
Nodup_MCC.base=MCC.base[!duplicated(MCC.base[,c("cm","cumI")]),]

#Sort by event dates#
sort_MCC.base=Nodup_MCC.base[order(Nodup_MCC.base[,"Time"]),]

#MCC is calculated by using the summation of cumulative incidences for all event of interest (first and recurrent)#
MCC=cumsum(sort_MCC.base[,"Deta"])

combine_MCC=cbind(sort_MCC.base,MCC)

#Only show the time points that have MCC change, and remove irrelevant information#
MCC.final=aggregate(combine_MCC[,"MCC"],list(combine_MCC[,"Time"]),max)###
colnames(MCC.final)=c("time","SumCIs")
} # for if(sum(cause==1)=0) condition.

return(MCC.final)
}
#-------------------------------------------------------------------------------------------------#


#-------------------------------------------------------------------------------------------------#
#---------------------- Sum of CIs with 95% confidence interval-----------------------------------#
SCI_95CI=function(id, time,cause,Tstart){
	
	MCC_org=SCI(id,time,cause,Tstart) ### run the original data to get MCC

	withdata=data.frame(id,time,cause,Tstart)
	uid=unique(withdata$id)
	for(boot in 1:1000){
		set.seed(boot+2016)
		sid=sample(uid,replace = TRUE)
		sid=data.frame(newid=seq(1,length(sid),1),id=sid)  ##use obs to be the new ids.
		bootdata=merge(withdata,sid,by.x="id",by.y="id")
		bootdata=data.frame(id=bootdata$newid,bootdata[, !names(bootdata) %in% c("id","newid")])
		#order the bootstrapped data  --no need, as the SCI function will order it inside the function ##
		#bootdata=bootdata[order(bootdata$id,bootdata$time,-bootdata$cause),]  
		out_boot=SCI(id=bootdata$id,time=bootdata$time,cause=bootdata$cause,Tstart=bootdata$Tstart)
		colnames(out_boot)=c("time",paste("MCC",boot,sep=""))
		MCC_org=merge(MCC_org,out_boot,by.x="time",by.y="time",all=T)
	}
	### For points not in a bootstrapped data, MCC is NA. Need to fill in by its previous MCC value ####.
	### If the first row is missing, the MCC=0
	MCC_org[1,][is.na(MCC_org[1,])]=0
	#FillDown: Fills in missing (NA) values with the previous non-missing value from "DataCombine" package.
#	for(boot in 1:1000){  ##this for loop takes a few minutes to run. Maybe use apply function to save time.
#		MCC_org=FillDown(MCC_org,paste("MCC",boot,sep=""))
#	} #"DataCombine" package not in R any more. Replace with the lines below. ==July 2025
	for(boot in 1:1000){  ##this for loop takes a few minutes to run. Maybe use apply function to save time.
		MCC_org[,paste("MCC",boot,sep="")]=na.locf(MCC_org[,paste("MCC",boot,sep="")])
	}
	### Take the 25th and 975th from each row to get 95% CI.
	sort_row=t(apply(MCC_org[,-c(1:2)],1,sort))
	MCC_95CI=data.frame(MCC_org[,1:2],low=sort_row[,25],up=sort_row[,975])
	#return(MCC_95CI)
	return(list("MCC_CI"=MCC_95CI,"MCCw1000vals"=MCC_org))
}


#-------------------------------------------------------------------------------------------------#
#-------------------------------------------------------------------------------------------------#



##--------------------------------------------------------------------------------------------------------------###
###### ---example data set that can be used for both left and right censoring, w or w/o 95% confidence intervals #######

dat0=read.csv("R:\\Biostatistics\\Biostatistics\\Common\\Qi\\webpage\\MCC\\dat for MCC14.csv")
dat0=dat0[order(dat0$subj),]

### If needed to svae time, we can only choose subset to test.
dat=dat0[1:600,]


dim(dat)
dat[1:5,]
id=dat$subj
time=dat$timedx  ###follow-up time from cohrot entry
cause=dat$status14
tstart=dat$astart  ##age at cohort entry
tstop=dat$timebirth  ##age at event/competing risk/censoring

########## right censonring only, no 95% CI -- takes a few seconds to run ##########
res1=MCC(id=id, time=time,cause=cause,type="SCI")  

######### left truncation and right censoring, no 95% CI -- takes a few seconds to run #######
res2=MCC(id=id, time=tstop,cause=cause,Tstart=tstart,type="SCI")  

 ##left truncation and right censoring, with 95% CI. --- takes ~20 minutes to run #######
res3=MCC(id=id, time=tstop,cause=cause,Tstart=tstart,type="SCI",c_interval=TRUE) 
names(res3)
### This is the MCC with 95% CI
MCC_CI=res3$MCC_CI
MCC_CI[1:10,]
### If the MCC and the 1000 bootstrpped values are needed, use
MCC1000=res3$MCCw1000vals
MCC1000[1:2,]


plot(c(0,70),c(0,8.5),type="n",xlab="Age",ylab="MCC")
 lines(MCC_CI$time,MCC_CI$SumCIs,type="s")
 lines(MCC_CI$time,MCC_CI$low,type="s",col=4)
 lines(MCC_CI$time,MCC_CI$up,type="s",col=4)



#####---------- example data: subsequent malignanct neoplamns (SMN), right censoring only --------------#####
##### It is a big datasets contains >14,000 individuals and hence bootstrap for 95% CI takes long to run ####

smn=read.csv("R:\\Biostatistics\\Biostatistics\\Common\\Qi\\webpage\\MCC\\smndata.csv")
smn[1:2,]
table(smn$SMNcat)
res_smn=MCC(id=smn$ccssid, time=smn$timesmn,cause=smn$SMNcat,type="SCI")   ## no 95% CI
res_smn95CI=MCC(id=smn$ccssid, time=smn$timesmn,cause=smn$SMNcat,type="SCI",c_interval=TRUE)   ##with 95% CI.1.5hrs
res_smn95CI=res_smn95CI$MCC_CI

plot(c(0,35),c(0,0.5),type="n",xlab="Years from DX",ylab="MCC")
 lines(res_smn95CI$time,res_smn95CI$SumCIs,type="s")
 lines(res_smn95CI$time,res_smn95CI$low,type="s",col=4)
 lines(res_smn95CI$time,res_smn95CI$up,type="s",col=4)

########## The "crprep" gave different results when R versions are different. Currently my tests showed that it ran as expected in R 3.5.X and 4.2.2 (but not 4.4.0). Users need to test the run whenever a new R version is used. Use this dataset to test: A correct version would give burden=1.3 at the last time point.
test_data=read.csv("R:\\Biostatistics\\Biostatistics\\Common\\Qi\\webpage\\MCC\\test_data.csv")
res_test=MCC(id=test_data$id, time=test_data$timebirth,cause=test_data$status15,Tstart=test_data$astart,type="SCI")